import numpy as np
import pandas as pd
import matplotlib.pylab as plt
import seaborn as sns
import math
from statsmodels.tsa.seasonal import seasonal_decompose
import statsmodels.api as sm
from pandas.plotting import lag_plot
from pandas.plotting import autocorrelation_plot
from statsmodels.graphics.tsaplots import plot_pacf
from statsmodels.tsa.stattools import adfuller
from itertools import cycle
from IPython.display import display
from math import ceil
color_pal = plt.rcParams['axes.prop_cycle'].by_key()['color']
color_cycle = cycle(plt.rcParams['axes.prop_cycle'].by_key()['color'])
# Preprocessing
from sklearn.preprocessing import MinMaxScaler
# Algorithms
from sklearn.cluster import KMeans
from sklearn.metrics import pairwise_distances_argmin_min
from sklearn.decomposition import PCA
from tslearn.barycenters import dtw_barycenter_averaging
print("setup complete")
sales = pd.read_csv('70prod_data.csv',\
index_col=0, parse_dates=True)
sales.head()
sales.describe().T
sales.plot(figsize=(20, 10),legend=False)
plt.show()
f, ax = plt.subplots(1, 1, figsize=(10, 10))
ax = sns.heatmap(sales.corr())
cg = sns.clustermap(sales.corr(), figsize=(18, 18), cbar_pos=(.1, .1, .03, .65))
cg.ax_row_dendrogram.set_visible(False)
cg.ax_col_dendrogram.set_visible(False)
plt.show()
On observe une forte correlation entre certains groupes de produits, ce qui suggère qu'on peut les classifier selon leurs comportements de ventes.
clustering¶
L'un des problèmes qu'on doit traiter avant de commencer le partitionnement (clustering) est l'échelle de la série. Sans normaliser les données, les séries qui se ressemblent seront vues si différentes les unes des autres et affecteront la précision du processus de clustering. Nous pouvons voir l'effet de la normalisation dans les images suivantes.
a = [[2],[7],[11],[14],[19],[23],[26]]
b = [[20000000],[40000000],[60000000],[80000000],[100000000],[120000000],[140000000]]
fig, axs = plt.subplots(1,3,figsize=(25,5))
axs[0].plot(a)
axs[0].set_title("Series 1")
axs[1].plot(b)
axs[1].set_title("Series 2")
axs[2].plot(a)
axs[2].plot(b)
axs[2].set_title("Series 1 & 2")
plt.figure(figsize=(25,5))
plt.plot(MinMaxScaler().fit_transform(a))
plt.plot(MinMaxScaler().fit_transform(b))
plt.title("Normalized Series 1 & Series 2")
plt.show()
mySeries = [sales[[col]] for col in sales.columns]
scaler = MinMaxScaler()
for i in range(len(mySeries)):
mySeries[i] = MinMaxScaler().fit_transform(mySeries[i])
mySeries[i]= mySeries[i].reshape(len(mySeries[i]))
Afin de regrouper nos séries avec des k-moyennes, les indices temporels des séries chronologiques seront considérés comme des caractéristiques différentes et seront les dimensions des points de données(les séries).
Puisque la longueur des série temporelles en pratique est souvent importante, un autre problème auquel on doit faire face est la malédiction de la dimensionnalité,ce terme a été inventé pour la première fois par Richard E. Bellman lors de l'examen des problèmes de programmation dynamique. Cela signifie essentiellement que lorsque la dimensionnalité des données augmente, la distance entre les points de données augmente également. Ainsi, ce changement de mesure de la distance affecte gravement les algorithmes basés sur la distance.
Pour résoudre ce problème, on va effectuer une analyse en composantes principales ACP avant le partitionnement.
pca = PCA(n_components=2)
pca_res = pca.fit_transform(mySeries)
Désormais avec moins de dimensions qu'avant, nous pouvons voir comment nos séries se répartissent en 2 dimensions.
plt.figure(figsize=(25,10))
plt.scatter(pca_res[:,0],pca_res[:,1], s=200)
plt.show()
print(pca_res[:10,:])
Kmeans: Le partitionnement en k-moyennes (ou k-means en anglais) est une méthode de partitionnement de données et un problème d'optimisation combinatoire. Étant donnés des points et un entier k, le problème est de diviser les points en k groupes, souvent appelés clusters, de façon à minimiser une certaine fonction. On considère la distance d'un point à la moyenne des points de son cluster (centroid) ; la fonction à minimiser est la somme des carrés de ces distances,la distance euclidienne est souvent utilisée.
Initialement on va éstimer le nombre de clusters requis par la racine du nombre des points (séries), on obtient 9 clusters.
cluster_count = math.ceil(math.sqrt(len(mySeries)))
cluster_count
et on va appliquer la méthode de Kmeans sur le résultat de l'ACP.
kmeans = KMeans(n_clusters=cluster_count,max_iter=5000)
labels = kmeans.fit_predict(pca_res)
Après l'entrainement du modèle, on a présenté les résultats. Pour chaque cluster, on a tracé les séries en gris, et afin de visualiser le mouvement et le comportement des séries appartenant au même cluster, on a pris la moyenne du cluster et tracé cette série moyenne en rouge.
som_x = som_y = math.ceil(math.sqrt(math.sqrt(len(mySeries))))
plot_count = math.ceil(math.sqrt(cluster_count))
fig, axs = plt.subplots(plot_count,plot_count,figsize=(25,25))
fig.suptitle('Clusters')
row_i=0
column_j=0
for label in set(labels):
cluster = []
for i in range(len(labels)):
if(labels[i]==label):
axs[row_i, column_j].plot(mySeries[i],c="gray",alpha=0.4)
cluster.append(mySeries[i])
if len(cluster) > 0:
#axs[row_i, column_j].plot(np.average(np.vstack(cluster),axis=0),c="red")
axs[row_i, column_j].plot(dtw_barycenter_averaging(np.vstack(cluster)),c="red")
axs[row_i, column_j].set_title("Cluster "+str(row_i*som_y+column_j))
column_j+=1
if column_j%plot_count == 0:
row_i+=1
column_j=0
plt.show()
Nous pouvons voir la distribution des séries chronologiques en clusters dans le graphique suivant.
les séries semblent uniformement distribuées aux clusters.
cluster_c = [len(labels[labels==i]) for i in range(cluster_count)]
cluster_n = ["cluster_"+str(i) for i in range(cluster_count)]
plt.figure(figsize=(15,5))
plt.title("Cluster Distribution for KMeans")
plt.bar(cluster_n,cluster_c)
plt.show()
Pour bénéficier de partitionnement qu'on a fait, il est nécessaire de lister pour chaque cluster, les séries qui lui appartiennent.
namesofMySeries=sales.columns
fancy_names_for_labels = [f"Cluster {label}" for label in labels]
groups = pd.DataFrame(zip(namesofMySeries,fancy_names_for_labels),columns=["Series","Cluster"]).\
sort_values(by="Cluster").set_index("Cluster")
for cluster in sorted(set(fancy_names_for_labels)):
print(cluster+':', \
', '.join(sorted(groups[groups.index == cluster]['Series'], \
key=lambda x: int(x.split('_')[1]))))
On a selectionné le produit le plus proche du centroide de chaque cluster, autrement dit le plus proche aux autres produits au même groupe.
closest, _ = pairwise_distances_argmin_min(kmeans.cluster_centers_, pca_res)
closest
Pour visualiser ce résultat, on a tracé l'ensemble des points ainsi que les centres des groupes (étoiles), le point le plus proche de chaque centre est présenté en gris.
np.unique(labels)
import matplotlib.cm as cm
plt.figure(figsize=(25,10))
plt.scatter(pca_res[:, 0], pca_res[:, 1], c=labels, s=300)
sc = plt.scatter(kmeans.cluster_centers_[:,0], kmeans.cluster_centers_[:,1], \
marker='*', c=np.unique(labels), edgecolors='black', s=500)
plt.scatter(pca_res[closest, 0], pca_res[closest, 1], c='grey', \
edgecolors=[sc.to_rgba(lab) for lab in np.unique(labels)], s=300)
for cluster_lab, cluster_x, cluster_y in \
zip(np.unique(labels), kmeans.cluster_centers_[:,0], kmeans.cluster_centers_[:,1]):
for pt_lab, pt_x, pt_y in zip(labels, pca_res[:, 0], pca_res[:, 1]):
if cluster_lab==pt_lab:
x, y = [cluster_x, pt_x], [cluster_y, pt_y]
plt.plot(x, y, c=sc.to_rgba(cluster_lab), alpha=0.5)
plt.show()
pour re vérifier ces résultats, on peut voir que les séries choisis sont centrées pour chaque groupe.
som_x = som_y = math.ceil(math.sqrt(math.sqrt(len(mySeries))))
plot_count = math.ceil(math.sqrt(cluster_count))
fig, axs = plt.subplots(plot_count,plot_count,figsize=(25,25))
fig.suptitle('Clusters')
row_i=0
column_j=0
for label in set(labels):
cluster = []
for i in range(len(labels)):
if(labels[i]==label):
if i not in closest :
axs[row_i, column_j].plot(mySeries[i], c="gray",alpha=0.4)
else :
axs[row_i, column_j].plot(mySeries[i],c="blue")
cluster.append(mySeries[i])
if len(cluster) > 0:
#axs[row_i, column_j].plot(np.average(np.vstack(cluster),axis=0),c="red")
axs[row_i, column_j].plot(dtw_barycenter_averaging(np.vstack(cluster)),c="red")
axs[row_i, column_j].set_title("Cluster "+str(row_i*som_y+column_j))
column_j+=1
if column_j%plot_count == 0:
row_i+=1
column_j=0
plt.show()
On a extrait les produits représentatifs des groupes et on a construit un nouveau jeu de données de séries temporelles sur lequel on va travailler par la suite.
sales_s = sales.iloc[:, sorted(closest)]
sales_s.head()
Puisque les produits choisis appartiennent à des groupes différents, on observe une grande dissemblance entre leurs caractéristiques statistiques.
display(sales_s.describe().T)
plt.figure(figsize=(10,5))
sns.boxplot(data=sales_s)
plt.show()
Ce qu'on peut confimer graphiquement :
sales_s.plot(subplots=True,layout=(3,3),figsize=(20,15), marker='.')
plt.show()
Moyenne mobile des séries¶
fig, axs = plt.subplots(5, 2, figsize=(20, 20))
axs = axs.flatten()
cols = [next(color_cycle), next(color_cycle), next(color_cycle)]
for i, item in enumerate(sales_s.columns):
sales_s[item].plot(title=item,
color=cols[0],
ax=axs[i],
marker='.', label="monthly")
sales_s[item].rolling(3).mean().plot(
color=cols[1],
ax=axs[i], label="quarter rolling mean", alpha=0.75)
sales_s[item].rolling(12).mean().plot(
color=cols[2],
ax=axs[i], linestyle='--', label="12months rolling mean")
axs[i].legend()
visualisation de la saisonnalité¶
fig, axes = plt.subplots(5, 2, figsize=(20, 20))
axes = axes.flatten()
for name, ax in zip(sales_s.columns, axes):
sns.boxplot(data=sales_s, x=sales_s.index.month, y=name, ax=ax)
Décomposition des séries¶
def plotseasonal(res, axes, name):
res.observed.plot(ax=axes[0], legend=False, title=name)
axes[0].set_ylabel('Observed')
res.trend.plot(ax=axes[1], legend=False)
axes[1].set_ylabel('Trend')
res.seasonal.plot(ax=axes[2], legend=False)
axes[2].set_ylabel('Seasonal')
res.resid.plot(ax=axes[3], legend=False, marker='o', linestyle='', markersize=5)
axes[3].set_ylabel('Residual')
plotseasonal.counter += 1
def grid_plotseasonal(df, ncol):
nrow=ceil(len(df.columns)/ncol)*4
fig, axes = plt.subplots(ncols=ncol, nrows=nrow, sharex=True, figsize=(15,1.5*nrow))
plotseasonal.counter, k = 0, 0
for i,col in enumerate(df.columns):
res = seasonal_decompose(df[col].dropna(), period=12, filt=None)
plotseasonal(res, axes[k:k+4, i%ncol], col)
if (i%ncol)==(ncol-1) : k+=4
plt.tight_layout()
plt.show()
grid_plotseasonal(sales_s, ncol=3)
La fonction season_decompose peut être utilisée pour l'analyse des portions de chaque composant de la série chronologique. Ceci est particulièrement utile pour déterminer l'absorption des résidus dans les données, sur la base des données décomposées. Le volume de cette absorption implique la prévisibilité de la série chronologique - plus les résidus sont élevés, moins la prévisibilité. Dans une certaine mesure, la proportion des résidus par rapport à la tendance et à la saisonnalité peut également être illustrée par les graphiques des moyennes mobiles et des écarts-types ci-dessus.
d = pd.DataFrame(0, index=sales_s.columns, columns=["RESMEAN","OBSMEAN","PERC"], dtype=float)
for col in sales_s.columns:
result = seasonal_decompose(sales_s[col], period=12, model='additive')
res, obs = result.resid, result.observed
d.loc[col][:2] = list(map(lambda x: np.mean(np.abs(x)), (res, obs[~np.isnan(res)])))
d.PERC = d.RESMEAN*100/d.OBSMEAN
d.round(2)
Un diagramme de décalage (lag plot) est un nuage de points pour une série chronologique et les mêmes données décalées. Avec un tel graphe, nous pouvons vérifier s'il existe une corrélation possible entre les ventes ce mois et le mois précédent.
fig, axes = plt.subplots(3, 3, figsize=(20, 15))
axes = axes.flatten()
for name, ax in zip(sales_s.columns, axes):
lag_plot(sales_s[name], ax=ax, label=name)
ax.legend()
Analyse d'autocorrélation¶
L'analyse d'autocorrélation illustre le potentiel de prédiction des données de séries chronologiques. Les graphiques d'autocorrélation résument graphiquement la force d'une relation avec une observation dans une série chronologique avec des observations à des pas de temps antérieurs. Le coefficient de Pearson est utilisé pour mesurer l'autocorrélation. Ainsi, l'analyse suivante n'est pertinente que pour les données avec une distribution gaussienne normale.
Un tracé de l'autocorrélation d'une série chronologique par décalage est appelé la fonction d'autocorrélation (ACF). Ce graphique est parfois appelé corrélogramme ou graphique d'autocorrélation. Le graphique montre la valeur de décalage le long de l'axe des x et la corrélation sur l'axe des y entre -1 et 1. Les intervalles de confiance sont dessinés sous forme de cône. Par défaut, il est défini sur un intervalle de confiance de 95%, ce qui suggère que les valeurs de corrélation en dehors de ce code sont très probablement une corrélation.
En général, la corrélation «partielle» entre deux variables est la quantité de corrélation entre elles qui n'est pas expliquée par leurs corrélations mutuelles avec un ensemble spécifié d'autres variables. Par exemple, si nous régressons une variable Y sur d'autres variables X1, X2 et X3, la corrélation partielle entre Y et X3 est la quantité de corrélation entre Y et X3 qui n'est pas expliquée par leurs corrélations communes avec X1 et X2.
fig, axes = plt.subplots(9, 2, figsize=(20, 35), sharex=False)
alpha=.05
for i, col in enumerate(sales_s.columns):
sm.graphics.tsa.plot_acf(sales[col].values.squeeze(), lags=40, ax=axes[i,0], title=f"{col} ACF", alpha=alpha)
sm.graphics.tsa.plot_pacf(sales[col].values.squeeze(), lags=15, ax=axes[i,1], title=f"{col} PACF", alpha=alpha)
# Import the check_stationarity function from previous lab
def stationarity_check(TS):
# Calculate rolling statistics
rolmean = TS.rolling(window = 8, center = False).mean()
rolstd = TS.rolling(window = 8, center = False).std()
# Perform the Dickey Fuller Test
dftest = adfuller(TS) # change the passengers column as required
#Plot rolling statistics:
fig = plt.figure(figsize=(12,6))
orig = plt.plot(TS, color='blue',label='Original')
mean = plt.plot(rolmean, color='red', label='Rolling Mean')
std = plt.plot(rolstd, color='black', label = 'Rolling Std')
plt.legend(loc='best')
plt.title('Rolling Mean & Standard Deviation')
plt.show(block=False)
# Print Dickey-Fuller test results
print ('Results of Dickey-Fuller Test:')
dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])
for key,value in dftest[4].items():
dfoutput['Critical Value (%s)'%key] = value
print(dfoutput)
return None
stationarity_check(sales_s.P_2)
différenciation d'ordre 1
fig, axs = plt.subplots(5, 2, figsize=(20, 20))
axs = axs.flatten()
for i, item in enumerate(sales_s.columns):
sales_s[item].diff().plot(title=f"{item} lag 1",
color=next(color_cycle),
ax=axs[i])
sales_s[item].diff().rolling(12).mean().plot(
color=cols[1],
ax=axs[i], label="rolling mean")
sales_s[item].diff().rolling(12).std().plot(
color=cols[2],
ax=axs[i], linestyle='--', label="rolling std")
axs[i].legend()
différenciation d'ordre 2, lag1 puis lag m
fig, axs = plt.subplots(5, 2, figsize=(20, 20))
axs = axs.flatten()
for i, item in enumerate(sales_s.columns):
sales_s[item].diff(1).diff(12).plot(title=f"{item} lag 1 lag 12",
color=next(color_cycle),
ax=axs[i])
sales_s[item].diff(1).diff(12).rolling(12).mean().plot(
color=cols[1],
ax=axs[i], label="rolling mean")
sales_s[item].diff(1).diff(12).rolling(12).std().plot(
color=cols[2],
ax=axs[i], linestyle='--', label="rolling std")
axs[i].legend()
